---
title: Evaluating the impact on F1 sterility of tsetse pupae Glossina morsitans morsitans irradiated following  short-term hypoxic conditioning 
author: "Mirieri et al.,"
date: "18-03-2025"
output: word_document
---

```{r }


#setwd("D:/2025/residual fertility_raw data_feb2025_2")

#load need library
library(datasets)
library(gcookbook)
library(ggplot2)
library(plyr)
library(dplyr)
library(lattice)
library(MASS)
library(rcompanion)
#library(survival)
library(ranger)
library(bbmle)
library(ggfortify)
library(knitr)
#library(coxme)
## Graphique "trellis"
library(lattice)
## Mod???¨les lin???©aires g???©n???©ralis???©s mixtes
library(lme4)
library(MuMIn) 
library(nlme)
## Graphics
library(tidyverse)
library(FSA)
library(stats)
library(RCA)
library(rmarkdown)
#library(broom)
library(sp)
library(ggpubr)
library(car)
library(ggthemes)
library(base)
```

```{r setup, include=FALSE, echo=FALSE}

require("knitr")
opts_knit$set(root.dir ="D:/2025/residual fertility_raw data_feb2025_2")
```

EXPERIMENT 1

```{r }

#####Duration on emergence and flight rates 


##emergence rates-combined hypoxia and irradiation
tab33 = read.csv("Fig1a.csv")
head(tab33)
Emerg_rate <- tab33$emerged / (tab33$unemerged+ tab33$emerged)
boxplot(Emerg_rate ~ tab33$Treatment_duration,xlab = "Treatment duration", ylab = "emergence rates")
boxplot(Emerg_rate ~ tab33$Irradiation_status,xlab = "Treatment duration", ylab = "emergence rates")


##relevel_relative significance

##1 hour
fm32 <- glmer(cbind(emerged, unemerged) ~ Treatments +(1|Replicate), family = binomial, data = tab33)
summary(fm32)

fm32a <- glmer(cbind(emerged, unemerged) ~ Irradiation_status +(1|Replicate), family = binomial, data = tab33)
summary(fm32a)
fm32b <- glmer(cbind(emerged, unemerged) ~ Irradiation_status_refIrrNorm +(1|Replicate), family = binomial, data = tab33)
summary(fm32b)

##No treatment
fm32 <- glmer(cbind(emerged, unemerged) ~ Treatment_duration +(1|Replicate), family = binomial, data = tab33)
summary(fm32)

##Irradiation in normoxia
fm33 <- glmer(cbind(emerged, unemerged) ~ Treatment_IN +(1|Replicate), family = binomial, data = tab33)
summary(fm33)
##chilling
fm34 <- glmer(cbind(emerged, unemerged) ~ Treatment_CH +(1|Replicate), family = binomial, data = tab33)
summary(fm34)
## no treatment
fm35 <- glmer(cbind(emerged, unemerged) ~ Treatment_NT +(1|Replicate), family = binomial, data = tab33)
summary(fm35)


### Figure 1a

Emerg_rate<- tab33$emerged / (tab33$unemerged+ tab33$emerged)
tab33$Treatment_duration<- as.factor(tab33$Treatment_duration)

            
Fig1a<-ggplot(tab33, aes(x=Treatment_duration,y=Emerg_rate, fill=Treatment_duration)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
              
Fig1a

tiff("Fig1a.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 600)
plot(Fig1a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Duration in hypoxia (hours)"))) + ylab(expression(bold("Emergence rate")))
dev.off()





####Flight rate_combined hypoxia and irradiation
tab1 = read.csv("Fig1b.csv")
head(tab1)
Flight_rate <- tab1$out / (tab1$in. + tab1$out)
boxplot(Flight_rate ~ tab1$Treatment_duration,xlab = "Duration in hypoxia", ylab = "Flight rates")
boxplot(Flight_rate ~ tab1$Irradiation_status,xlab = "Duration in hypoxia", ylab = "Flight rates")

### Figure1b

Flight_rate<- Flight_rate <- tab1$out / (tab1$in. + tab1$out)
tab1$Treatment_duration<- as.factor(tab1$Treatment_duration)

            
Fig1b<-ggplot(tab1, aes(x=Treatment_duration,y=Flight_rate, fill=Treatment_duration)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
              
              
Fig1b

tiff("Fig1b.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 600)
plot(Fig1b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Duration in hypoxia (hours)"))) + ylab(expression(bold("Flight propensity")))
dev.off()




##relevel_relative significance
##1hour
fm1a <- glmer(cbind(out,in.) ~ Treatment +(1|Replicate), family = binomial, data = tab1)
summary(fm1a)

##no treatment
fm1a <- glmer(cbind(out,in.) ~ Treatment_duration +(1|Replicate), family = binomial, data = tab1)
summary(fm1a)

##normoxia
fm1 <- glmer(cbind(out,in.) ~ Treatment_IN +(1|Replicate), family = binomial, data = tab1)
summary(fm1)

##chilling
fm2 <- glmer(cbind(out,in.) ~ Treatment_CH +(1|Replicate), family = binomial, data = tab1)
summary(fm2)
###No treatment
fm3 <- glmer(cbind(out,in.) ~ Treatment_NT +(1|Replicate), family = binomial, data = tab1)
summary(fm3)

fm4 <- glmer(cbind(out,in.) ~ Irradiation_status +(1|Replicate), family = binomial, data = tab1)
summary(fm4)

fm4a <- glmer(cbind(out,in.) ~ Irradiation_status_refIrrNorm +(1|Replicate), family = binomial, data = tab1)
summary(fm4a)



```



EXPERIMENT 2
```{r }

####Dose response curve for mating 
tab1 = read.csv("Fig2.csv")

tab1$Irradiation_dose=as.factor(tab1$Irradiation_dose)

head(tab1)
mate <- tab1$pairs_formed/ (tab1$unformed_pairs+ tab1$pairs_formed)
boxplot(mate ~ tab1$Irradiation_doses,xlab = "Treatment (dose)", ylab = "Mating ability")
boxplot(mate ~ tab1$Irradiation_status,xlab = "Irradiation conditions", ylab = "Mating ability")



#wrapped according to irradiation conditions 
ggplot(tab1, aes(x=Irradiation_doses,y=mate, fill=Irradiation_status)) +
  geom_boxplot(alpha=0.0) +
  geom_jitter(width=0.1,alpha=0.2) +
  xlab("Pupal_age")+ 
  facet_wrap(~Irradiation_status,ncol = 4) +
  theme(axis.text.x = element_text(angle = 0, hjust = 1))
#theme(axis.text.x = element_text(angle = 45, hjust = 1))

Fig5<-ggplot(tab1,aes(x=Irradiation_status,y=mate, fill=Irradiation_dose)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)



fm1a <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_dose +(1|Replicate), family = binomial, data = tab1)
summary(fm1a)


###relevel
tab1$Irradiation_dose <- relevel(tab1$Irradiation_dose, ref = "110")
fm4 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_dose +(1|Replicate), family = binomial, data = tab1)
summary(fm4)


fm1b <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_status +(1|Replicate), family = binomial, data = tab1)
summary(fm1b)


fm1c <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_dose + Irradiation_status+(1|Replicate), family = binomial, data = tab1)
summary(fm1c)


tab1$Irradiation_dose <- relevel(tab1$Irradiation_dose, ref = "0")
fm4 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_dose+ Irradiation_status +(1|Replicate), family = binomial, data = tab1)
summary(fm4)

fm1d <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_status * Irradiation_status+(1|Replicate), family = binomial, data = tab1)
summary(fm1d)



#########
#boxplot(mate ~ tab1$Irradiation_doses, xlab ="Irradiation dose", ylab = "Mating ability")
#boxplot(mate ~ tab1$Irradiation._status, xlab ="Irradiation status", ylab = "Mating ability")


### mating ability within hypoxia (pairwise comparisons and significance)
hyp<-subset(tab1, Irradiation_status=="Hypoxia")
summary(hyp)
mate <- hyp$pairs_formed/ (hyp$unformed_pairs+ hyp$pairs_formed)
boxplot(mate ~ hyp$Irradiation_dose, xlab ="Irradiation dose (Gy)", ylab = "Mating ability (Hypoxia)")

mod9 <- glm(mate~Irradiation_doses+(1|Replicate), data = hyp)
summary(mod9)


###relevel
hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "110")
fm4 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_dose +(1|Replicate), family = binomial, data = hyp)
summary(fm4)

###relevel
hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "150")
fm4 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_dose +(1|Replicate), family = binomial, data = hyp)
summary(fm4)



##Mating ability within normoxia (pairwise comparisons and significance)

norm<-subset(tab1, Irradiation_status=="Normoxia")
mate <- norm$pairs_formed/ (norm$unformed_pairs+ norm$pairs_formed)
boxplot(mate ~ norm$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Mating ability (Normoxia)")


mod9 <- glm(mate~Irradiation_doses+(1|Replicate), data =norm)
summary(mod9)


###relevel
norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "110")
fm4 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_dose +(1|Replicate), family = binomial, data = norm)
summary(fm4)

###relevel
norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "150")
fm4 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Irradiation_dose +(1|Replicate), family = binomial, data = norm)
summary(fm4)




###Mating ability -dose response curve wrapped according to irradiation conditions 

#Figure 2
mate <- tab1$pairs_formed/ (tab1$unformed_pairs+ tab1$pairs_formed)

Fig2<-ggplot(tab1,aes(x=Irradiation_status2,y=mate, fill=Irradiation_dose2)) +
  geom_boxplot(position = position_dodge(preserve = 'single')) 
Fig2


tiff("Fig2.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig2+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold(" Mating ability ")))
dev.off()


```


EXPERIMENT 2

```{r }

###Dose response curve on survival of males and females (pairwise boxplots and significance)


data_4<- read.csv("Fig3a_3b.csv")
summary(data_4)
head(data_4)
data_4$Irradiation_dose=as.factor(data_4$Irradiation_dose)


Survival_21days<-data_4$surv21

boxplot(Survival_21days ~ data_4$NT, ylab = "Survival rate after 21 days")

boxplot(Survival_21days ~ data_4$Irradiation_status, ylab = "Survival rate after 21 days")
boxplot(Survival_21days ~ data_4$sex._surv, ylab = "Survival rate after 21 days")


fm4_1 <- glmer(cbind(surv21,n - surv21) ~ NT +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

fm4_1 <- glmer(cbind(surv21,n - surv21) ~ Irradiation_status +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

fm4_1 <- glmer(cbind(surv21,n - surv21) ~ Irradiation_status2 +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

fm4_1 <- glmer(cbind(surv21,n - surv21) ~ sex._surv +(1|rep), family = binomial, data = data_4)
summary(fm4_1)


#*#*#*#Figure 3a


## survival of male and females wrapped according to irradiation status 
data_4$Irradiation_doses=as.factor(data_4$Irradiation_doses)

Fig3a<-ggplot(data_4,aes(x=Irradiation_status2,y=Survival_21days, fill=Irradiation_doses)) +
  geom_boxplot(position = position_dodge(preserve = 'single')) 
Fig3a

tiff("Fig3a.tiff", width = 6, height = 4, units = 'in', res = 300)
plot(Fig3a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Survival for 21 days")))
dev.off()



#*#*#*#Figure 3b

## survival of flies wrapped according to sex (Male and Female) 
Fig3b<-ggplot(data_4,aes(x=sex._surv2,y=Survival_21days, fill=Irradiation_doses)) +
  geom_boxplot() 
Fig3b

tiff("Fig3b.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig3b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Sex of flies"))) + ylab(expression (bold("Survival for 21 days ")))
dev.off()


####survival of male and female within the hypoxic conditions 

shyp<-subset(data_4 , Irradiation_status=="Hypoxia")
head(shyp)
Survival_21days<-shyp$surv21

boxplot(Survival_21days ~ shyp$Irradiation_doses, ylab = "Survival rate after 21 days(hypoxia)")
boxplot(Survival_21days ~ shyp$sex._surv, ylab = "Survival rate after 21 days(hypoxia)")


fm5_1 <- glmer(cbind(surv21,n - surv21) ~ NT +(1|rep), family = binomial, data = shyp)
summary(fm5_1)

Survival_21days<-shyp$surv21

fm5_1 <- glmer(cbind(surv21,n - surv21) ~ sex._surv +(1|rep), family = binomial, data = shyp)
summary(fm5_1)

Survival_21days<-shyp$surv21



#### survival of males and females within the normoxic conditions 
snorm<-subset(data_4 , Irradiation_status=="Normoxia")
head(snorm)

Survival_21days<-snorm$surv21

boxplot(Survival_21days ~ snorm$Irradiation_doses, ylab = "Survival rate after 21 days(Normoxia)")
boxplot(Survival_21days ~ snorm$sex._surv, ylab = "Survival rate after 21 days(Normoxia)")


###sig
fm5_1 <- glmer(cbind(surv21,n - surv21) ~ treatment +(1|rep), family = binomial, data = snorm)
summary(fm5_1)
Survival_21days<-snorm$surv21


fm5_1 <- glmer(cbind(surv21,n - surv21) ~ sex._surv +(1|rep), family = binomial, data = snorm)
summary(fm5_1)

Survival_21days<-snorm$surv21



```

EXPERIMENT 2

```{r }

### Pupae Per Initial Female (PPIF (fecundity))- dose response curve (pairwise comparisons of boxplots and significance ) 


Prod = read.csv("Fig4a.csv")



Prod$Irradiation_dose=as.factor(Prod$Irradiation_dose)



head(Prod)
summary(Prod)

boxplot(Prod$PPIF ~ Prod$Irradiation_doses, xlab ="Irradiation dose(Gy)", ylab = "PPIF_dose response")
boxplot(Prod$PPIF ~ Prod$Irradiation_status, xlab ="Irradiation status", ylab = "PPIF_dose response")



#fm1 <- glm(PPIF~ Irradiation_dose + Irradiation_status, data = Prodres)
#fm2 <- glm(PPIF~Irradiation_dose * Irradiation_status, data = Prodres)
#fm3 <- glm(PPIF~ Irradiation_status, data = Prod)
#fm4<- glm(PPIF~ Irradiation_dose, data = Prodres)


#fm1 <- glm(PPIF~ relevelNT_Irradiation_doses  + Irradiation_status, data = Prodres)
#fm2 <- glm(PPIF~relevelNT_Irradiation_doses  * Irradiation_status, data = Prodres)
fm3 <- glm(PPIF~ Irradiation_status, data = Prod)
fm3b <- glm(PPIF~ Irradiation_status2, data = Prod)
fm4 <- glm(PPIF~ relevelNT_Irradiation_doses , data = Prod)


#summary(fm1)
#summary(fm2)
summary(fm3)
summary(fm3b)
summary(fm4)



### Pupae production within hypoxic conditions 

hyp<-subset(Prod, Irradiation_status=="Hypoxia")
summary(hyp)
boxplot(hyp$PPIF ~ hyp$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "PPIF_dose response(Hypoxia)")


mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

###relevel
hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "0")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "70")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "90")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "110")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "130")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)


hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "150")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "170")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "190")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)




###Pupae production within hypoxic conditions normoxic conditions 

norm<-subset(Prod, Irradiation_status=="Normoxia")
summary(norm)
boxplot(norm$PPIF ~ norm$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "PPIF_dose response (Normoxia)")


mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

###relevel
norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "0")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "70")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "90")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "110")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "130")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "150")
mod9 <- glm(PPIF~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

#*#*#Figure 4a


### All in one graph- PPIF (fecundity)
Fig4a<-ggplot(Prod,aes(x=Irradiation_status2,y=PPIF, fill=Irradiation_dose)) +
  geom_boxplot(position = position_dodge(preserve = 'single'))
Fig4a

tiff("Fig4a.tiff", width =6, height = 4, units = 'in', res = 600)
plot(Fig4a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("PPIF")))
dev.off()




##Residual fertility_ dose response curve

res = read.csv("Fig4b.csv")

res$Irradiation_dose=as.factor(res$Irradiation_dose)


boxplot(res$Residual_fertlity ~ res$relevelNT_Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Residual fertlity_dose response")
boxplot(res$Residual_fertlity ~ res$Irradiation_status, xlab ="Irradiation status", ylab = "Residual fertlity_dose response")



fm1 <- glm(Residual_fertlity~ Irradiation_dose + Irradiation_status, data = res)
fm2 <- glm(Residual_fertlity~Irradiation_dose * Irradiation_status, data = res)
fm8 <- glm(Residual_fertlity~ Irradiation_status, data = res)
fm9 <- glm(Residual_fertlity~ Irradiation_dose, data = res)

summary(fm1)
summary(fm2)
summary(fm8)
summary(fm9)


####residual fertility within hypoxic conditions 
hyp<-subset(res, Irradiation_status=="Hypoxia")
summary(hyp)
boxplot(hyp$Residual_fertlity ~ hyp$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Residual fertlity_dose response (Hypoxia)")

mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

###relevel
hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "0")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "70")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "90")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "110")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "130")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)


hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "150")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "170")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

hyp$Irradiation_dose <- relevel(hyp$Irradiation_dose, ref = "190")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = hyp)
summary(mod9)

### Residual fertility within normoxic conditions 
norm<-subset(res, Irradiation_status=="Normoxia")
summary(norm)
boxplot(norm$Residual_fertlity ~ norm$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Residual fertlity_dose response (Normoxia)")


mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

###relevel
norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "0")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "70")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "90")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "110")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "130")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

norm$Irradiation_dose <- relevel(norm$Irradiation_dose, ref = "150")
mod9 <- glm(Residual_fertlity~Irradiation_dose+(1|Replicate), data = norm)
summary(mod9)

#*##Figure 4b


###All in one graph -residual fertility
Fig4b<-ggplot(res,aes(x=Irradiation_status2,y=Residual_fertlity, fill=Irradiation_dose)) +
  geom_boxplot() 
Fig4b

tiff("Fig4b.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig4b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Residual fertility")))
dev.off()


```

EXPERIMENT 3

```{r }

###PUPUAE PRODUCTION (fecundity)

## Pupae production-  pairwise boxplot comparisons and significance

data_7<- read.csv("Fig5a_5bsig.csv")
head(data_7)
summary(data_7)
data_7=na.omit(data_7)

#**#
boxplot(data_7$ Pupae_production~ data_7$ Irradiation_status,xlab = " Irradiation status", ylab = "Pupae_production")
boxplot(data_7$ Pupae_production~ data_7$ Irradiation_dose,xlab = " Irradiation dose", ylab = "  Pupae_production")
boxplot(data_7$ Pupae_production~ data_7$ Parental_sex_fromtreatment,xlab = " Parental_sex_fromtreatment ", ylab = "Pupae_production")
boxplot(data_7$ Pupae_production~ data_7$ Generation,xlab = " Generation", ylab = "  Pupae_production")



fm1<- glm(Pupae_production ~Irradiation_dose+(1|Rep), data = data_7)
summary(fm1)
fm2 <- glm(Pupae_production~Irradiation_status+(1|Rep), data = data_7)
summary(fm2)

##relevel
fm3 <- glm(Pupae_production~Irradiation_status2+(1|Rep), data = data_7)
summary(fm3)

fm4<- glm(Pupae_production ~Parental_sex_fromtreatment+(1|Rep), data = data_7)
summary(fm4)
fm5 <- glm(Pupae_production~Generation_res+(1|Rep), data = data_7)
summary(fm5)

####models 
mod1<- glm(Pupae_production ~Irradiation_dose+Irradiation_status+Parental_sex_fromtreatment+(1|Rep), data = data_7)
mod2<- glm(Pupae_production ~Irradiation_dose+Irradiation_status+Generation_res+(1|Rep), data = data_7)
mod3<- glm(Pupae_production ~Irradiation_dose+Irradiation_status+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_7)
mod4<- glm(Pupae_production~Irradiation_dose+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_7)
mod5<- glm(Pupae_production ~Irradiation_status+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_7)
mod6<- glm(Pupae_production ~Irradiation_status+Generation_res+(1|Rep), data = data_7)
mod5<- glm(Pupae_production~Irradiation_status+Parental_sex_fromtreatment+(1|Rep), data = data_7)
mod7<- glm(Pupae_production~Irradiation_dose+Irradiation_status++(1|Rep), data = data_7)
mod8<- glm(Pupae_production~Irradiation_dose+(1|Rep), data = data_7)
mod9 <- glm(Pupae_production~Irradiation_status+(1|Rep), data = data_7)
mod10 <- glm(Pupae_production~Irradiation_status2+(1|Rep), data = data_7)
mod11<- glm(Pupae_production ~Parental_sex_fromtreatment+(1|Rep), data = data_7)
mod12<- glm(Pupae_production~Generation_res+(1|Rep), data = data_7)

AICc(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9,mod10,mod11,mod12)

summary(mod3)




### all treatments in one graph
data_7<- read.csv("Fig5a.csv")

head(data_7)
summary(data_7)
data_7=na.omit(data_7)

#*#*##Figure 5a

###Pupae production-all treatments in one graph
boxplot(data_7$Pupae_production~ data_7$ Treatments,xlab = " Treatments", ylab = "Pupae_production")

Fig5a<-ggplot(data_7,aes(x=Treatments2,y=Pupae_production, fill=Treatments2)) +
  geom_boxplot() 
Fig5a

tiff("Fig5a.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig5a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) +  theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("All treatments"))) + ylab(expression (bold("Pupae production")))
dev.off()




data_7<- read.csv("Fig5b.csv")

###Pupae production- all treatment wrapped according to irradiation conditions 
boxplot(data_7$Pupae_production~ data_7$ Treatments,xlab = " Treatments", ylab = " Pupae_production")


#*#*##Figure 5b


Fig5b<-ggplot(data_7,aes(x=Irradiation_status2,y=Pupae_production, fill=Treatments2)) +
  geom_boxplot()
Fig5b

tiff("Fig5b.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig5b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) +  theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Pupae production")))
dev.off()

fm8<- glm(Pupae_production~Irradiation_status+(1|Replicate), data = data_7)
summary(fm8)



### Pupae Per Initial Female (PPIF)


###all treatments -pairwise comparisons (boxplots and significance)

data_7<- read.csv("Fig5c_5dsig.csv")
head(data_7)
summary(data_7)
data_7=na.omit(data_7)


boxplot(data_7$ PPIF~ data_7$ Irradiation_status,xlab = " Irradiation status", ylab = "PPIF")
boxplot(data_7$ PPIF~ data_7$ Irradiation_dose,xlab = " Irradiation dose", ylab = "  PPIF")
boxplot(data_7$ PPIF~ data_7$ Parental_sex_fromtreatment,xlab = " Parental_sex_fromtreatment ", ylab = "PPIF")
boxplot(data_7$ PPIF~ data_7$ Generation,xlab = " Generation", ylab = "PPIF")



fm1<- glm(PPIF ~Irradiation_dose+(1|Rep), data = data_7)
summary(fm1)
fm2 <- glm(PPIF~Irradiation_status+(1|Rep), data = data_7)
summary(fm2)

##relevel
fm3 <- glm(PPIF~Irradiation_status2+(1|Rep), data = data_7)
summary(fm3)

fm4<- glm(PPIF ~Parental_sex_fromtreatment+(1|Rep), data = data_7)
summary(fm4)
fm5 <- glm(PPIF~Generation_res+(1|Rep), data = data_7)
summary(fm5)

####models 
mod1<- glm(PPIF ~Irradiation_dose+Irradiation_status+Parental_sex_fromtreatment+(1|Rep), data = data_7)
mod2<- glm(PPIF ~Irradiation_dose+Irradiation_status+Generation_res+(1|Rep), data = data_7)
mod3<- glm(PPIF ~Irradiation_dose+Irradiation_status+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_7)
mod4<- glm(PPIF~Irradiation_dose+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_7)
mod5<- glm(PPIF ~Irradiation_status+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_7)
mod6<- glm(PPIF ~Irradiation_status+Generation_res+(1|Rep), data = data_7)
mod5<- glm(PPIF~Irradiation_status+Parental_sex_fromtreatment+(1|Rep), data = data_7)
mod7<- glm(PPIF~Irradiation_dose+Irradiation_status++(1|Rep), data = data_7)
mod8<- glm(PPIF~Irradiation_dose+(1|Rep), data = data_7)
mod9 <- glm(PPIF~Irradiation_status+(1|Rep), data = data_7)
mod10 <- glm(PPIF~Irradiation_status2+(1|Rep), data = data_7)
mod11<- glm(PPIF ~Parental_sex_fromtreatment+(1|Rep), data = data_7)
mod12<- glm(PPIF~Generation_res+(1|Rep), data = data_7)

AICc(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9,mod10,mod11,mod12)

summary(mod3)



### all treatments on one graph-PPIF
data_8<- read.csv("Fig5c.csv")
head(data_8)
summary(data_8)
data_8=na.omit(data_8)

###PPIF 
boxplot(data_8$PPIF~ data_8$ Treatments1,xlab = " Treatments", ylab = "PPIF")

#*#*## Figure 5c 
Fig5c<-ggplot(data_8,aes(x=Treatments1,y=PPIF, fill=Treatments1)) +
  geom_boxplot() 
Fig5c

tiff("Fig5c.tiff", width = 9, height = 4, units = 'in', res = 600)
plot(Fig5c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("All treatments"))) + ylab(expression (bold("PPIF")))
dev.off()



#### All treatments -PPIF (boxplots and significance)
data_8<- read.csv("Fig5d.csv")
head(data_8)
summary(data_8)
data_8=na.omit(data_8)


###each irradiation condition
hypx1<-subset(data_8, Irradiation_status=="Hypoxia")

summary(hypx1)
head(hypx1)


boxplot(hypx1$PPIF~ hypx1$Treatments,xlab = " Treatments", ylab = " PPIF")



norm1<-subset(data_8, Irradiation_status=="Normoxia ")
summary(norm1)
head(norm1)


boxplot(norm1$PPIF~ norm1$Treatments,xlab = " Treatments", ylab = " PPIF")


#*#*## Figure 5d

###### PPIF-all treatments wrapped according irradiation conditions

Fig5d<-ggplot(data_8,aes(x=Irradiation_status2,y=PPIF, fill=Treatment1)) +
  geom_boxplot()
Fig5d

tiff("Fig5d.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig5d+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("PPIF")))
dev.off()

fm8<- glm(PPIF ~Irradiation_status+(1|Rep), data = data_8)
summary(fm8)

```



EXPERIMENT 3
```{r}

###Emergence rates for F1 AND F2

data_6<- read.csv("Fig6a_sig.csv")
head(data_6)
summary(data_6)

#Figure 5

###emergence rates
data_6$rate <- data_6$emerged / (data_6$unemerged+ data_6$emerged)

boxplot(data_6$rate ~ data_6$ Irradiation_status,xlab = " Irradiation status", ylab = "Emergence rates")
boxplot(data_6$rate ~ data_6$ Irradiation_dose,xlab = " Irradiation_dose", ylab = "Emergence rates")
boxplot(data_6$rate ~ data_6$ Sex_treatment_parental,xlab = " Sex_treatment_parental", ylab = "Emergence rates")
boxplot(data_6$rate ~ data_6$ Emerging.pupae.generation,xlab = " Emerging.pupae.generation", ylab = "Emergence rates")


#boxplot(data_6$rate ~ data_6$ Parental_generation,xlab = " Parental_generation", ylab = "Emergence rates")



fm1a <- glmer(cbind(emerged, unemerged) ~ Irradiation_status +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm1a <- glmer(cbind(emerged, unemerged) ~ Irradiation_dose +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm1a <- glmer(cbind(emerged, unemerged) ~ Sex_treatment_parental +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)


fm1a <- glmer(cbind(emerged, unemerged) ~ Emerging.pupae.generation +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)


###emergence rates all treatments in one graph


data_6b<- read.csv("Fig6a.csv")
head(data_6b)
summary(data_6b)
###emergence rates all in one  graph 
data_6b$rate <- data_6b$emerged / (data_6b$unemerged+ data_6b$emerged)

boxplot(data_6b$rate ~ data_6b$ Treatments,xlab = " Treatments", ylab = "Emergence rates")

Fig6a<-ggplot(data_6b,aes(x=Treatments2,y=data_6b$rate, fill=Treatments2)) +
  geom_boxplot() 
Fig6a


tiff("Fig6a.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig6a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold(" All treatments"))) + ylab(expression (bold("Emergence rate")))
dev.off()



###sex ratios all in one graph separated into normoxia and hypoxia


### male emergence rates all in one graph

data_6b<- read.csv("Fig7a_7c.csv")
head(data_6b)
summary(data_6b)
###emergence rates all in one  graph 



data_6b$male<- data_6b$Male_emergence / (data_6b$Female_emergence+ data_6b$Male_emergence)

boxplot(data_6b$male~ data_6b$Treatments,xlab = " Treatments", ylab = " Male Emergence rates")


Fig7a<-ggplot(data_6b,aes(x=Treatments2,y=data_6b$male, fill=Treatments2)) +
  geom_boxplot() 
Fig7a

tiff("Fig7a.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig7a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold(" All treatments "))) + ylab(expression (bold("Male emergence rate")))
dev.off()

### female emergence rates all in one graph
data_6b$female<- data_6b$Female_emergence / (data_6b$Male_emergence+ data_6b$Female_emergence)

boxplot(data_6b$female~ data_6b$ Treatments,xlab = " Treatments", ylab = " Female Emergence rates")

Fig7c<-ggplot(data_6b,aes(x=Treatments2,y=data_6b$female, fill=Treatments2)) +
  geom_boxplot() 
Fig7c

tiff("Fig7c.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig7c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("All treatments"))) + ylab(expression (bold("Female emergence rate")))
dev.off()


### emergence, all treatments wrapped according to irradiation conditions (hypoxia or normoxia)

data_6b<- read.csv("Fig6b.csv")
head(data_6b)
summary(data_6b)

###emergence rates all in one graph wrapped according to irradiation conditions
exp3emerge <- data_6b$emerged / (data_6b$unemerged+ data_6b$emerged)

Fig6b<-ggplot(data_6b,aes(x=Irradiation_status2,y=exp3emerge, fill=Treatment1)) +
  geom_boxplot() 
Fig6b

tiff("Fig6b.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig6b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Emergence rate")))
dev.off()






###### Emergence sex ratios as influenced by different treatments individually

###Male emergence (sex of parental generation from treatment )-FO had only males

data_6$male<- data_6$Male_emergence / (data_6$Female_emergence+ data_6$Male_emergence)

boxplot(data_6$male~ data_6$ Irradiation_status,xlab = " Irradiation status", ylab = " Male Emergence rates")
boxplot(data_6$male~ data_6$ Irradiation_dose ,xlab = " Irradiation status", ylab = " Male Emergence rates")
boxplot(data_6$male~ data_6$ Emerging.pupae.generation ,xlab = " Emerging.pupae.generation", ylab = " Male Emergence rates")



fm1a <- glmer(cbind(Male_emergence, Female_emergence) ~ Irradiation_status +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm1a <- glmer(cbind(Male_emergence, Female_emergence) ~ Irradiation_dose +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm1a <- glmer(cbind(Male_emergence, Female_emergence) ~ Emerging.pupae.generation  +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)


###Female emergence (sex of parental generation )

data_6$female<- data_6$Female_emergence / (data_6$Male_emergence+ data_6$Female_emergence)

boxplot(data_6$female~ data_6$ Irradiation_status,xlab = " Irradiation status", ylab = " Female Emergence rates")
boxplot(data_6$female~ data_6$ Irradiation_dose,xlab = " Irradiation dose", ylab = " Female Emergence rates")
boxplot(data_6$female~ data_6$ Emerging.pupae.generation,xlab = " Emerging.pupae.generation", ylab = " Female Emergence rates")



fm1 <- glmer(cbind(Female_emergence, Male_emergence) ~ Irradiation_status +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm2 <- glmer(cbind(Female_emergence, Male_emergence) ~ Irradiation_dose +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm3 <- glmer(cbind(Female_emergence, Male_emergence) ~ Emerging.pupae.generation +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)



###sex ratios- individual boxplots and significance (separated by irradiation conditions)

##Normoxia

norm<-subset(data_6, Irradiation_status=="Normoxia")

norm$female<- norm$Female_emergence / (norm$Male_emergence+ norm$Female_emergence)
norm$male<- norm$Male_emergence / (norm$Female_emergence+ norm$Male_emergence)


boxplot(norm$female ~norm$Emerging.pupae.generation, xlab =" Emergence pupae generation", ylab = "Emergence (normoxia-female)")

boxplot(norm$female ~norm$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Emergence(normoxia-female)")

boxplot(norm$male ~norm$Emerging.pupae.generation, xlab =" Emergence pupae generation", ylab = "Emergence (normoxia-male)")

boxplot(norm$male ~norm$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Emergence(normoxia-male)")


fm1nf<- glm(norm$female ~Emerging.pupae.generation+(1|Replicate), data = norm)
summary(fm1nf)

fm2nf<- glm(norm$female ~Irradiation_dose+(1|Replicate), data = norm)
summary(fm2nf)

fm1nm<- glm(norm$male ~Emerging.pupae.generation+(1|Replicate), data = norm)
summary(fm1nm)

fm2nm<- glm(norm$male ~Irradiation_dose+(1|Replicate), data = norm)
summary(fm2nm)




######### Hypoxia
hypx<-subset(data_6, Irradiation_status=="Hypoxia")

hypx$female<- hypx$Female_emergence / (hypx$Male_emergence+ hypx$Female_emergence)
hypx$male<- hypx$Male_emergence / (hypx$Female_emergence+ hypx$Male_emergence)


boxplot(hypx$female ~hypx$Emerging.pupae.generation, xlab =" Emergence pupae generation", ylab = "Emergence (Hypoxia-female)")

boxplot(hypx$female ~hypx$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Emergence(Hypoxia- female)")

boxplot(hypx$male ~hypx$Emerging.pupae.generation, xlab =" Emergence pupae generation", ylab = "Emergence (Hypoxia-male)")

boxplot(hypx$male ~hypx$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Emergence(Hypoxia-male)")


fm1hf<- glm(hypx$female ~Emerging.pupae.generation+(1|Replicate), data = hypx)
summary(fm1hf)

fm2hf <- glm(hypx$female ~Irradiation_dose+(1|Replicate), data = hypx)
summary(fm2hf)
fm1hm<- glm(hypx$male ~Emerging.pupae.generation+(1|Replicate), data = hypx)
summary(fm1hm)

fm2hm <- glm(hypx$male ~Irradiation_dose+(1|Replicate), data = hypx)
summary(fm2hm)




###sex ratios all in one graph separated according to irradiation conditions
data_6b<- read.csv("Fig7d_7b.csv")
head(data_6b)
summary(data_6b)



#females (normoxia/hypoxia )
data_6b$female<- data_6b$Female_emergence / (data_6b$Male_emergence+ data_6b$Female_emergence)

Fig7d<-ggplot(data_6b,aes(x=Irradiation_status2,y=data_6b$female, fill=Treatments)) +
  geom_boxplot() 
Fig7d

tiff("Fig7d.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig7d+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Female emergence rate ")))
dev.off()



#males (normoxia/hypoxia )
data_6b$male<- data_6b$Male_emergence / (data_6b$Female_emergence+ data_6b$Male_emergence)

Fig7b<-ggplot(data_6b,aes(x=Irradiation_status2,y=data_6b$male, fill=Treatments)) +
  geom_boxplot() 
Fig7b

tiff("Fig7b.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig7b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Male emergence rate")))
dev.off()


```


EXPERIMENT 3

```{r }

###Emergence rates for F1 AND F2-pairwise comparison boxplots and significance

data_6<- read.csv("Fig6a_sig.csv")
head(data_6)
summary(data_6)


###emergence rates
data_6$rate <- data_6$emerged / (data_6$unemerged+ data_6$emerged)

boxplot(data_6$rate ~ data_6$ Irradiation_status,xlab = " Irradiation status", ylab = "Emergence rates")
boxplot(data_6$rate ~ data_6$ Irradiation_dose,xlab = " Irradiation_dose", ylab = "Emergence rates")
boxplot(data_6$rate ~ data_6$ Sex_treatment_parental,xlab = " Sex_treatment_parental", ylab = "Emergence rates")
boxplot(data_6$rate ~ data_6$ Emerging.pupae.generation,xlab = " Emerging.pupae.generation", ylab = "Emergence rates")


#boxplot(data_6$rate ~ data_6$ Parental_generation,xlab = " Parental_generation", ylab = "Emergence rates")



fm1a <- glmer(cbind(emerged, unemerged) ~ Irradiation_status +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm1a <- glmer(cbind(emerged, unemerged) ~ Irradiation_dose +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm1a <- glmer(cbind(emerged, unemerged) ~ Sex_treatment_parental +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)


fm1a <- glmer(cbind(emerged, unemerged) ~ Emerging.pupae.generation +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)



###emergence rates all treatments in one graph


data_6b<- read.csv("Fig6a.csv")
head(data_6b)
summary(data_6b)

data_6b$rate <- data_6b$emerged / (data_6b$unemerged+ data_6b$emerged)

boxplot(data_6b$rate ~ data_6b$ Treatments,xlab = " Treatments", ylab = "Emergence rates")



#*#*##Figure 6a

Fig6a<-ggplot(data_6b,aes(x=Treatments2,y=data_6b$rate, fill=Treatments2)) +
  geom_boxplot() 
Fig6a

tiff("Fig6a.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig6a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold(" All treatments"))) + ylab(expression (bold("Emergence rate")))
dev.off()


### emergence, all in one graph wrapped according to irradiation conditions (hypoxia or normoxia)

data_6b<- read.csv("Fig6b.csv")
head(data_6b)
summary(data_6b)

###emergence rates all in one graph
exp3emerge <- data_6b$emerged / (data_6b$unemerged+ data_6b$emerged)

#*#*## Figure 6b
Fig6b<-ggplot(data_6b,aes(x=Irradiation_status2,y=exp3emerge, fill=Treatment1)) +
  geom_boxplot() 
Fig6b

tiff("Fig6b.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig6b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Emergence rate")))
dev.off()






### male emergence rates all in one graph

data_6b<- read.csv("Fig7a_7c.csv")
head(data_6b)
summary(data_6b)
###emergence rates all in one  graph 



data_6b$male<- data_6b$Male_emergence / (data_6b$Female_emergence+ data_6b$Male_emergence)

boxplot(data_6b$male~ data_6b$Treatments,xlab = " Treatments", ylab = " Male Emergence rates")



#*#**# Figure 7a
Fig7a<-ggplot(data_6b,aes(x=Treatments2,y=data_6b$male, fill=Treatments2)) +
  geom_boxplot() 
Fig7a

tiff("Fig7a.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig7a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold(" All treatments "))) + ylab(expression (bold("Male emergence rate")))
dev.off()

### female emergence rates all in one graph
data_6b$female<- data_6b$Female_emergence / (data_6b$Male_emergence+ data_6b$Female_emergence)

boxplot(data_6b$female~ data_6b$ Treatments,xlab = " Treatments", ylab = " Female Emergence rates")


#*#*## Figure 7c


Fig7c<-ggplot(data_6b,aes(x=Treatments2,y=data_6b$female, fill=Treatments2)) +
  geom_boxplot() 
Fig7c

tiff("Fig7c.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig7c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("All treatments"))) + ylab(expression (bold("Female emergence rate")))
dev.off()



###### Emergence sex - pairwise comparison boxplots and significance

###Male emergence (sex of parental generation from treatment )-FO had only males irradiated

data_6$male<- data_6$Male_emergence / (data_6$Female_emergence+ data_6$Male_emergence)

boxplot(data_6$male~ data_6$ Irradiation_status,xlab = " Irradiation status", ylab = " Male Emergence rates")
boxplot(data_6$male~ data_6$ Irradiation_dose ,xlab = " Irradiation status", ylab = " Male Emergence rates")
boxplot(data_6$male~ data_6$ Emerging.pupae.generation ,xlab = " Emerging.pupae.generation", ylab = " Male Emergence rates")



fm1a <- glmer(cbind(Male_emergence, Female_emergence) ~ Irradiation_status +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm1a <- glmer(cbind(Male_emergence, Female_emergence) ~ Irradiation_dose +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm1a <- glmer(cbind(Male_emergence, Female_emergence) ~ Emerging.pupae.generation  +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)


###Female emergence (sex of parental generation )

data_6$female<- data_6$Female_emergence / (data_6$Male_emergence+ data_6$Female_emergence)

boxplot(data_6$female~ data_6$ Irradiation_status,xlab = " Irradiation status", ylab = " Female Emergence rates")
boxplot(data_6$female~ data_6$ Irradiation_dose,xlab = " Irradiation dose", ylab = " Female Emergence rates")
boxplot(data_6$female~ data_6$ Emerging.pupae.generation,xlab = " Emerging.pupae.generation", ylab = " Female Emergence rates")



fm1 <- glmer(cbind(Female_emergence, Male_emergence) ~ Irradiation_status +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm2 <- glmer(cbind(Female_emergence, Male_emergence) ~ Irradiation_dose +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)

fm3 <- glmer(cbind(Female_emergence, Male_emergence) ~ Emerging.pupae.generation +(1|Replicate), family = binomial, data = data_6)
summary(fm1a)



###sex ratios- individual boxplots and significance according to irradiation conditions 

##Normoxic conditions 

norm<-subset(data_6, Irradiation_status=="Normoxia")

norm$female<- norm$Female_emergence / (norm$Male_emergence+ norm$Female_emergence)
norm$male<- norm$Male_emergence / (norm$Female_emergence+ norm$Male_emergence)


boxplot(norm$female ~norm$Emerging.pupae.generation, xlab =" Emergence pupae generation", ylab = "Emergence (normoxia-female)")

boxplot(norm$female ~norm$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Emergence(normoxia-female)")

boxplot(norm$male ~norm$Emerging.pupae.generation, xlab =" Emergence pupae generation", ylab = "Emergence (normoxia-male)")

boxplot(norm$male ~norm$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Emergence(normoxia-male)")


fm1nf<- glm(norm$female ~Emerging.pupae.generation+(1|Replicate), data = norm)
summary(fm1nf)

fm2nf<- glm(norm$female ~Irradiation_dose+(1|Replicate), data = norm)
summary(fm2nf)

fm1nm<- glm(norm$male ~Emerging.pupae.generation+(1|Replicate), data = norm)
summary(fm1nm)

fm2nm<- glm(norm$male ~Irradiation_dose+(1|Replicate), data = norm)
summary(fm2nm)




######### Hypoxic condition 
hypx<-subset(data_6, Irradiation_status=="Hypoxia")

hypx$female<- hypx$Female_emergence / (hypx$Male_emergence+ hypx$Female_emergence)
hypx$male<- hypx$Male_emergence / (hypx$Female_emergence+ hypx$Male_emergence)


boxplot(hypx$female ~hypx$Emerging.pupae.generation, xlab =" Emergence pupae generation", ylab = "Emergence (Hypoxia-female)")

boxplot(hypx$female ~hypx$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Emergence(Hypoxia- female)")

boxplot(hypx$male ~hypx$Emerging.pupae.generation, xlab =" Emergence pupae generation", ylab = "Emergence (Hypoxia-male)")

boxplot(hypx$male ~hypx$Irradiation_dose, xlab ="Irradiation dose(Gy)", ylab = "Emergence(Hypoxia-male)")


fm1hf<- glm(hypx$female ~Emerging.pupae.generation+(1|Replicate), data = hypx)
summary(fm1hf)

fm2hf <- glm(hypx$female ~Irradiation_dose+(1|Replicate), data = hypx)
summary(fm2hf)
fm1hm<- glm(hypx$male ~Emerging.pupae.generation+(1|Replicate), data = hypx)
summary(fm1hm)

fm2hm <- glm(hypx$male ~Irradiation_dose+(1|Replicate), data = hypx)
summary(fm2hm)



###sex ratios all in one graph separated according to irradiation conditions
data_6b<- read.csv("Fig7d_7b.csv")
head(data_6b)
summary(data_6b)

#*#*## Figure 7d
#females (normoxia/hypoxia )
data_6b$female<- data_6b$Female_emergence / (data_6b$Male_emergence+ data_6b$Female_emergence)

Fig7d<-ggplot(data_6b,aes(x=Irradiation_status2,y=data_6b$female, fill=Treatments)) +
  geom_boxplot() 
Fig7d

tiff("Fig7d.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig7d+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Female emergence rate ")))
dev.off()


#*#*## Figure 7b

#males (normoxia/hypoxia )
data_6b$male<- data_6b$Male_emergence / (data_6b$Female_emergence+ data_6b$Male_emergence)

Fig7b<-ggplot(data_6b,aes(x=Irradiation_status2,y=data_6b$male, fill=Treatments)) +
  geom_boxplot() 
Fig7b

tiff("Fig7b.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig7b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Male emergence rate")))
dev.off()


```



EXPERIMENT 3

```{r }

###  Residual fertility of F0 and F1 combined- Pairwise comparisons of all treatments (boxplots and significance) 

data_8<- read.csv("Fig8sig.csv")

data_8
head(data_8)
summary(data_8)


boxplot(data_8$Residual_fertility~ data_8$ Irradiation_status,xlab = " Irradiation status", ylab = " Residual fertility")
boxplot(data_8$Residual_fertility~ data_8$ Irradiation_dose,xlab = " Irradiation dose", ylab = " Residual fertility")
boxplot(data_8$Residual_fertility~ data_8$ Parental_sex_fromtreatment,xlab = " Parental_sex_fromtreatment ", ylab = " Residual fertility")
boxplot(data_8$Residual_fertility~ data_8$ Generation_res,xlab = " Generation", ylab = " Residual fertility")
#boxplot(data_8$Residual_fertility~ data_8$ parental_sex_gen,xlab = " parental sex and generation", ylab = " Residual fertility")


######
fm1<- glm(Residual_fertility ~Irradiation_dose+(1|Rep), data = data_8)
summary(fm1)
fm2 <- glm(Residual_fertility ~Irradiation_status+(1|Rep), data = data_8)
summary(fm2)
fm3<- glm(Residual_fertility ~Parental_sex_fromtreatment+(1|Rep), data = data_8)
summary(fm3)
fm4 <- glm(Residual_fertility ~Generation_res+(1|Rep), data = data_8)
summary(fm4)
#fm8 <- glm(Residual_fertility ~parental_sex_gen+(1|Rep), data = data_8)
#summary(fm8)

####models 
mod1<- glm(Residual_fertility ~Irradiation_dose+Irradiation_status+Parental_sex_fromtreatment+(1|Rep), data = data_8)
mod2<- glm(Residual_fertility ~Irradiation_dose+Irradiation_status+Generation_res+(1|Rep), data = data_8)
mod3<- glm(Residual_fertility ~Irradiation_dose+Irradiation_status+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_8)
mod4<- glm(Residual_fertility ~Irradiation_dose+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_8)
mod5<- glm(Residual_fertility ~Irradiation_status+Parental_sex_fromtreatment+Generation_res+(1|Rep), data = data_8)
mod6<- glm(Residual_fertility~Irradiation_status+Generation_res+(1|Rep), data = data_8)
mod5<- glm(Residual_fertility ~Irradiation_status+Parental_sex_fromtreatment+(1|Rep), data = data_8)
mod7<- glm(Residual_fertility ~Irradiation_dose+Irradiation_status++(1|Rep), data = data_8)
mod8<- glm(Residual_fertility ~Irradiation_dose+(1|Rep), data = data_8)
mod9 <- glm(Residual_fertility~Irradiation_status+(1|Rep), data = data_8)
mod10 <- glm(Residual_fertility~Irradiation_status+(1|Rep), data = data_8)
mod11<- glm(Residual_fertility ~Parental_sex_fromtreatment+(1|Rep), data = data_8)
mod12<- glm(Residual_fertility~Generation_res+(1|Rep), data = data_8)

AICc(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9,mod10,mod11,mod12)
##best is mod7
#anova(mod1, mod7)

summary(mod3)


### Residual fertility-All treatments in one graph

data_8<- read.csv("Fig8a.csv")

head(data_8)
summary(data_8)

Fig8a<-ggplot(data_8,aes(x=Treatment,y=Residual_fertility, fill=Treatment)) +
  geom_boxplot()
Fig8a

tiff("Fig8a.tiff", width = 7, height = 4, units = 'in', res = 600)
plot(Fig8a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) +theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("All treatments"))) + ylab(expression (bold("Residual fertility")))
dev.off()



######Generations as affected by irradiation dose and conditions

data_8<- read.csv("Fig8sig.csv")


######## Residual fertility of F0-Pairwise comparisons (boxplots and significance)

genf0<-subset(data_8, Generation_res=="F0")
genf0
###boxplot(genf0$Residual_fertility~ genf0$ Parental_sex_fromtreatment,xlab = " parental sex and generation", ylab = " Residual fertility")
boxplot(genf0$Residual_fertility~ genf0$ Irradiation_dose,xlab = " Irradiation_dose", ylab = " Residual fertility")
boxplot(genf0$Residual_fertility~ genf0$ Irradiation_status,xlab = " Irradiation_status", ylab = " Residual fertility")

## not there##fm8a<- glm(Residual_fertility ~Parental_sex_fromtreatment+(1|Rep), data = genf0)
###summary(fm8a)
fm8b<- glm(Residual_fertility ~Irradiation_dose+(1|Rep), data = genf0)
summary(fm8b)
fm8c <- glm(Residual_fertility ~Irradiation_status+(1|Rep), data =genf0)
summary(fm8c)

#fm8 <- glm(Residual_fertility ~parental_sex_gen+(1|Rep), data = genf1)
#summary(fm8)



###### Residual feritility of F1-Pairwise comparisons (boxplots and significance)


genf1<-subset(data_8, Generation_res=="F1")
genf1
boxplot(genf1$Residual_fertility~ genf1$ Parental_sex_fromtreatment,xlab = " parental sex and generation", ylab = " Residual fertility")
boxplot(genf1$Residual_fertility~ genf1$ Irradiation_dose,xlab = " Irradiation_dose", ylab = " Residual fertility")
boxplot(genf1$Residual_fertility~ genf1$ Irradiation_status,xlab = " Irradiation_status", ylab = " Residual fertility")


fm8d <- glm(Residual_fertility ~Parental_sex_fromtreatment+(1|Rep), data = genf1)
summary(fm8d)

fm8e<- glm(Residual_fertility ~Irradiation_dose+(1|Rep), data = genf1)
summary(fm8e)
fm8f <- glm(Residual_fertility ~Irradiation_status+(1|Rep), data = genf1)
summary(fm8f)

#fm8 <- glm(Residual_fertility ~parental_sex_gen+(1|Rep), data = data_8)
#summary(fm8)




### Residual fertility according to irradiation conditions 

###Normoxic conditions 
norm1<-subset(data_8, Irradiation_status=="Normoxia ")
head(norm1)
str(norm1)

##
boxplot(norm1$Residual_fertility~ norm1$ Irradiation_dose,xlab = " Irradiation dose", ylab = " Residual fertility")
boxplot(norm1$Residual_fertility~ norm1$ Parental_sex_fromtreatment,xlab = " Parental_sex_fromtreatment ", ylab = " Residual fertility")
boxplot(norm1$Residual_fertility~ norm1$ Generation,xlab = " Generation", ylab = " Residual fertility")



####
fm8<- glm(Residual_fertility ~Irradiation_dose+(1|Rep), data = norm1)
summary(fm8)
fm8<- glm(Residual_fertility ~Parental_sex_fromtreatment+(1|Rep), data = norm1)
summary(fm8)
fm8 <- glm(Residual_fertility ~Generation_res+(1|Rep), data = norm1)
summary(fm8)

#### Hypoxic conditions 

hyp1<-subset(data_8, Irradiation_status=="Hypoxia")
head(hyp1)
str(hyp1)

##
boxplot(hyp1$Residual_fertility~ hyp1$ Irradiation_dose,xlab = " Irradiation dose", ylab = " Residual fertility")
boxplot(hyp1$Residual_fertility~ hyp1$ Parental_sex_fromtreatment,xlab = " Parental_sex_fromtreatment ", ylab = " Residual fertility")
boxplot(hyp1$Residual_fertility~ hyp1$ Generation,xlab = " Generation", ylab = " Residual fertility")

###
fm8<- glm(Residual_fertility ~Irradiation_dose+(1|Rep), data = hyp1)
summary(fm8)
fm8<- glm(Residual_fertility ~Parental_sex_fromtreatment+(1|Rep), data = hyp1)
summary(fm8)
fm8 <- glm(Residual_fertility ~Generation_res+(1|Rep), data = hyp1)
summary(fm8)




###Residual fertility- all treatments (wrapped according to irradiation conditions) 

data_8<- read.csv("Fig8b.csv")
head(data_8)
summary(data_8)
data_8=na.omit(data_8)


##for figure

boxplot(data_8$Residual_fertility~ data_8$ Irradiation_status,xlab = " Treatments ", ylab = " Residual fertility(hypoxia)")

boxplot(data_8$Residual_fertility~ data_8$ Irradiation_status,xlab = " Treatments", ylab = " Residual fertility")

Fig8b<-ggplot(data_8,aes(x=Irradiation_status,y=Residual_fertility, fill=Treatments1)) +
  geom_boxplot()
Fig8b

tiff("Fig8b.tiff", width = 6, height = 4, units = 'in', res = 600)
plot(Fig8b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_text(colour = "black")) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation conditions"))) + ylab(expression (bold("Residual fertility")))
dev.off()

fm8<- glm(PPIF ~Irradiation_status+(1|Rep), data = data_8)
summary(fm8)



```
